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ABSTRACT 

Over the past thirty years, numerical methods and simulation tools for incompressible flows have been 
advanced as a subset of the computational fluid dynamics (CFD) discipline. Although incompressible flows are 
encountered in many areas of engineering, simulation of compressible flow has been the major driver for 
developing computational algorithms and tools. This is probably due to the rather stringent requirements for 
predicting aerodynamic performance characteristics of flight vehicles, while flow devices involving low-speed 
or incompressible flow could be reasonably well designed without resorting to accurate numerical simulations. 
As flow devices are required to be more sophisticated and highly efficient, CFD tools become increasingly 
important in fluid engineering for incompressible and low-speed flow. This paper reviews some of the successes 
made possible by advances in computational technologies during the same period, and discusses some of the 
current challenges faced in computing incompressible flows. 

L INTRODUCTION 

Almost all, if not all, fluid is compressible. Incompressible flow is an approximation of flow where flow speed 
is insignificant eveiywhere compared to the speed of sound of the medium. If incompressible flow is defined in 
this way, the majority of the fluid and associated flow we encounter in our daily lives belong to the 
incompressible category. For example, most of the flow associated with air and water (such as flow related to 
automobiles, ships, submarines, the water supply through pipes and channels, hydraulic turbines, pumps, low- 
speed airplanes, and trout swimming in mountain streams) and flow of biofluid (such as blood) are all in the 
incompressible flow domain. One of the earliest mathematical models of incompressible flow is the famous 
equation by Bernoulli, who in 1730 developed the model equation while investigating blood flow (see an 
interesting article by Quarteroni [1]). It is not surprising that scientists have been investigating incompressible 
flow analytically, experimentally, and computationally ever since. 

Computational study of incompressible flow problems in both basic research and engineering applications has 
been performed for several decades. Numerical solutions for such basic fluid dynamics problems as: flow past a 
circular cylinder; flow through channels, ducts and pipes; and flow over a backward facing step were presented 
as early as the 1930s (e.g., Thom [2] for a circular cylinder). In the computational fluid dynamics (CFD) 
community, especially in aerospace, CFD is synonymous with computational aerodynamics. Computational 
incompressible fluid dynamics was treated as a special case of aerodynamics, that is, the low-speed limit of 
compressible fluid dynamics. This is probably due to requirements to produce highly accurate results for 
predicting aerodynamic performance characteristics of flight vehicles, while flow devices involving low-speed 
flow can be reasonably well designed without resorting to accurate numerical simulations (e.g., hydraulic 
turbines for hydroelectric power plants were designed without the CFD approach). As flow devices become 
increasingly compact and efficient — pushing the conventional operating envelope — requirements on 
incompressible CFD tools have become more demanding, just as aerodynamic performance prediction tools 
require quantitative prediction capability. This trend is reflected in the development of various incompressible 
flow solution methods and tools, especially, in conjunction with high-fidelity computations using high-end 
computing facilities. 

Mathematically, the incompressible flow formulation poses unique issues not present in compressible equations 
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because of the incompressibility requirement. Physically, information travels at infinite speed in an 
incompressible medium, which imposes stringent requirements on computational algorithms for satisfying 
incompressibility as well as difficulties in designing downstream boundaiy conditions. Differences in various 
methods of solving the incompressible flow equations originate from differences in strategies of satisfying 
incompressibility. 

Generally, two different approaches are used in solving viscous incompressible flow equations, i.e., 
incompressible Navier-Stokes equations. The first is based on satisfying the incompressibility directly. If 
primitive variables, i.e., pressure and velocity, are chosen, pressure is used as a mapping parameter to satisfy 
incompressibility. This class of methods is generally known as the “pressure-based” method. Instead of pressure 
and velocity, derived quantities such as stream function- vorticity and vorticity- velocity can be used, resulting in 
different sets of governing equations. For general three-dimensional applications, however, the primitive 
variable formulation poses the least difficulty in geometry modeling and in setting the boundary conditions. The 
second approach is based on compressible flow formulation where momentum and continuity equations are 
coupled through the use of density. Incompressibility is recovered as a limiting case of this formulation. This 
class of methods is known as the “density-based” method. The artificial compressibility method is typical of 
this approach, which is a special case of preconditioned compressible flow methods. 

During the past several decades, many review articles and books on CFD have discussed incompressible flow 
methods. For more comprehensive reviews of computational methods for incompressible flow in general, 
readers are referred to these materials, i.e.: Roach [3]; Peyret and Taylor [4]; Hirsch [5]; Kwak [6]; Gunzburger 
and Nicolades [7]; Hafez and Oshima [8]; Gresho and Sani [9]; and Hafez [10]. These books and articles 
provide fairly extensive formulations, numerical methods, and solutions to ftmdamental fluid dynamic 
problems. This paper gives a short summary of methods from a historical perspective. After thirty-some years 
of numerical simulation of incompressible flows, industrial problems involving complex systems are now 
solved routinely. This paper focuses on engineering applications to review what has been accomplished to date 
and what challenges remain to be resolved. There are a vast number of cases where incompressible CFD 
methods have made significant impacts. Those presented here represent limited samples to illustrate the level of 
complexities researchers encounter in fluid engineering. 

IL EVOLUTION OF SOLUTION METHODS 

For general applications, the approach based on primitive variables is the most flexible choice. Other 
formulations using derived variables, such as stream function and vorticity, eliminate the need for pressure 
boundary conditions. However, this comes at the expense of adding a new requirement for derived variables, 
such as vorticity boundary conditions. All these formulations produce comparable results for two-dimensional 
problems, however, extension to three dimensions is not straightforward for stream function formulation. For 
vorticity- velocity formulation in three dimensions, three Poisson-type vorticity equations are introduced that 
make the formulation very expensive to solve. In this paper, the examples of three-dimensional applications 
have primarily been obtained using primitive variable formulation. Therefore, this section briefly reviews the 
evolution of the two most popular approaches using the primitive variables. Other formulations and examples 
are found in the literature [3-10]. 

2.1 Formulation 

Three-dimensional incompressible flow with constant density is governed by the following Navier-Stokes 
equations: 


0 

dxi 

du~ du-'Uj dp dr tj 

dt dXj dx t dXj 


( 2 . 1 ) 

( 2 . 2 ) 


2 


Where t is the time, x { the Cartesian coordinates, Ui the corresponding velocity components, p the pressure, and 
Tg the viscous stress tensor. Here, all variables are non-dimensionalized by a reference velocity and length 

scale. For turbulent flow, T t - should include turbulent stresses. 

The major difference between the incompressible and the compressible Navier-Stokes formulations is in the 
continuity equation. The incompressible formulation can be viewed as a singular limit of the compressible one. 
Satisfying the mass conservation equation, therefore, is the primary issue in solving the above set of equations. 
Physically, incompressible flow is characterized by elliptic behavior of the pressure waves, the speed of which 
in a truly incompressible flow is infinite. The pressure field is desired as a part of the solution; however, the 
pressure condition is needed on the boundary for numerical computation. One can utilize other formulations 
using derived quantities like vorticity to eliminate the pressure from the boundary conditions; however, 
boundary conditions for the derived variables have to be imposed. In realistic 3-D problems, these derived 
quantities are difficult to define or impractical to use. The primitive -variable formulation, namely, using 
pressure and velocities as dependent variables, then becomes very convenient and flexible in 3-D applications. 
However, in this formulation, mass conservation and its relation to pressure must be properly handled while 
achieving computational efficiency. Although various techniques have been developed in the past, none have 
proven to be universally better than the others. 

2.2 Methods Based on Pressure Iteration 


In 1965, Harlow and Welch [11] published the first primitive variable method using a derived Poisson equation 
for pressure. In this method, called the marker-and-cell (MAC) method, the pressure is used as a mapping 
parameter to satisfy the continuity equation. By taking the divergence of the momentum equation, the Poisson 
equation for pressure is obtained: 
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The usual computational procedure involves choosing the pressure field at the current time step such that 
continuity is satisfied at the next time step. The original MAC method is based on a staggered arrangement on a 
2-D Cartesian grid. The staggered grid conserves mass, momentum, and kinetic energy in a natural way and 
avoids odd-even point decoupling of the pressure encountered in a regular grid [12]. However, these properties 
become unclear when applied to generalized curvilinear coordinates. Even though the original method used an 
explicit Euler solver, various time-advancing schemes can be implemented here. Ever since its introduction, 
numerous variations of the MAC method have been devised and successful computations have been made. 

The major drawback of this method is the large amount of computing time required for solving the Poisson 
equation for pressure. When the physical problem requires a very small time step, the penalty paid for an 
iterative solution procedure for the pressure may be tolerable. But the method as a whole is slow and the 
pressure boundary condition is difficult to specify. One important aspect of the numerical solution of the 
Poisson equation for pressure is tied to the spatial differencing of the second derivatives. To satisfy the mass 
conservation in grid space, the difference form of the second derivative in the Poisson equation has to be 
constructed consistently with the discretized momentum equation (see Kwak [6] for more detail). 

For a steady-state solution, the correct pressure field is desired only when the solution is converged, so the 
iteration procedure for the pressure can be simplified such that it requires only a few iterations at each time step. 
The best-known method using this approach is the Semi-Implicit Method for Pressure-Linked Equations 
(SIMPLE) [13-14]. The unique feature of this method is the simple way of estimating the velocity correction 
and thus the pressure correction. This feature simplifies the computation but introduces empiricism into the 
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method. Despite its empiricism, the SIMPLE method has been used successfully for many computations, and 
has been the basis for many pressure-based methods developed since (e.g., [15] and many others). It is not 
however, the intention of this paper to fully review the evolution of this method. Many successful applications 
and tools developed to date are the real testimonial of the value of this pressure-based approach, in general. 

Fractional-step procedure for pressure-based method 

The MAC method can be viewed as a special case of the projection method (i.e., Chorin, [16]) or, from the 
operator point of a view, one variation of the fractional-step method [17-19]). Here, operator splitting can be 
accomplished in several ways by treating the momentum equation as a combination of convection, pressure, and 
viscous terms. The common application of this method is done in two steps: The first step is to solve for an 
auxiliary velocity field using the momentum equation, in which the pressure- gradient term can be computed 
from the pressure in the previous time step or can be excluded entirely. In the second step, the pressure is 
computed, which can map the auxiliary velocity onto a divergence-free velocity field. This procedure was first 
applied in staggered Cartesian grids by Harlow and Welch [11], and later extended to regular grids (e.g., Dwyer 
et al. [20] and many others). 

One particular aspect of the fractional-step method requiring special care is the intermediate boundary 
condition. Orszag et al. [21] discussed this extensively. Rosenfeld et al. [22] devised a generalized scheme 
where physical boundary conditions can be used at intermediate steps. As with other pressure-based methods, 
the efficiency of the fractional-step method depends on the Poisson solver. A multi-grid acceleration, which is 
physically consistent with the elliptic field, is one possible avenue to enhance the computational efficiency. In 
light of succftssfi.il computations in Cartesian coordinates using staggered grid, it became of interest to develop a 
method for general three-dimensional coordinates. However, extension of the staggered arrangement to general 
computational cells is not straightforward and involves features related to geometric conserv ation. A general 
numerical method based on a staggered grid was developed at NASA Ames Research Center (among others), 
initially by Rosenfeld et al. [22]. Later, Kiris and Kwak [23] developed a more robust implicit procedure for 
“not so smooth grids” using a finite-volume framework on a staggered grid. A staggered grid has favorable 
properties in Cartesian coordinates, such as coupling odd-even points. However, it is debatable whether a 
staggered grid has clear advantages over a regular grid in generalized curvilinear coordinates. In practice, this 
is because flux balancing in general curved cells requires interpolation that makes the differences between a 
cell-centered and staggered arrangement insignificant. Further discussion related to this approach is found m 
recent collection of papers by Hafez [10] (for example, in papers by Loner et al. [24] and more). 

2.3 Artificial Compressibility Method 

Advances in CFD for aerospace have been made in conjunction with compressible flow computations. 
Therefore, it is of significant interest to be able to use some of these compressible flow algorithms. One of the 
early attempts to do this was the artificial compressibility method by Chorin (1967) [25]. In this formulation, 
the continuity equation is modified by adding a time-derivative of the pressure term resulting in: 


1 dp + du t 
dt dx t 


(2.4) 


where fi is an artificial compressibility or a pseudo-compressibility parameter. These two terms, artificial 


compressibility and pseudo-compressibility, are used interchangeably in the literature. Together with the 
unsteady momentum equations, this forms a hyperbolic-parabolic type of time-dependent system of equations. 
Thus, the implicit schemes developed for compressible flows can be implemented. Note that t no longer 
represents a true physical time in this formulation. 


Physically, this means that waves of finite speed are introduced into the incompressible flow 7 field as a medium 
to distribute the pressure. For a truly incompressible flow, the wave speed is infinite, whereas the speed of 
propagation of these pseudo-waves depends on the magnitude of the artificial compressibility. Additionally, in 
a true incompressible flow, the pressure field is affected instantaneously by any disturbances in the flow field, 
but with artificial compressibility, there is a time lag between the flow disturbance and its effect on the pressure 
field. Ideally, the value of the artificial compressibility is chosen to be as high as the particular choice of 
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algorithm allows, so that the incompressibility is recovered quickly. This must be done without decreasing the 
accuracy and the stability property of the numerical method implemented. On the other hand, if the artificial 
compressibility is chosen such that these waves travel too slowly, then the variation of the pressure field 
accompanying these waves will be very slow. This will interfere with the timely development of the viscous 
boundary layer. In viscous flows, the behavior of the boundary layer is very sensitive to the streamwise pressure 
gradient, especially when the boundary layer is separated. If separation is present, a pressure wave traveling 
with finite speed will cause a change in the local pressure gradient, which will affect the location of the flow 
separation. This change in separated flow will feed back to the pressure field, possibly preventing convergence 
to a steady state. Especially for internal flow, the viscous effect is important for the entire flow field and the 
interaction between the pseudo pressure-waves and the viscous flow field becomes veiy important. 

Artificial compressibility relaxes the strict requirement to satisfy mass conservation in each step. However, to 
utilize this convenient feature, it is essential to understand the nature of artificial compressibility both physically 
and mathematically. Chang and Kwak (1984) [26] reported details of artificial compressibility, and suggested 
some useful guidelines for choosing the artificial compressibility parameter. 

For obtaining steady state solutions, the governing equations are iterated in pseudo-time until the solution 
converges to steady state, at which time the original incompressible Navier-Stokes equations are satisfied. 
Time-dependent calculation of incompressible flows are especially time consuming due to the elliptic nature of 
the governing equations. Physically, this means that any local change in the flow has to be propagated to the 
entire flow field. Numerically, this means that in each time step, the pressure field has to go through one 
complete steady-state iteration cycle; for example, by a Poisson-solver-type pressure iteration or an artificial 
compressibility iteration method. In transient flow, the physical time step has to be small and, consequently, the 
change in the flow field may be small. In this situation, the number of iterations in each time step for obtaining 
a divergence-free flow field may not be as high as for regular steady-state computations. However, time- 
accurate computations are in general extremely time-consuming. Therefore, it is particularly desirable to 
develop computationally efficient methods by implementing a fast algorithm and by utilizing computer 
characteristics such as parallel processing. 

Various applications that evolved from the artificial compressibility concept have been reported for obtaining 
steady-state solutions (e.g., Kwak et al. [27], Chang et al. [28], Choi and Merkle [29]). To obtain time- 
dependent solutions using this method, an iterative procedure can be applied in each physical time step such that 
the continuity equation is satisfied. Merkle and Athavale [30] and Rogers et al. [31] reported successful 
computations using this pseudo-time iteration approach. 

Pre-conditioning 

In incompressible flow, the density is constant; this is the source of difficulty in solving the compressible 
Navier-Stokes equation as the flow approaches the incompressible limit. Adding a non-physical pressure term 
[25] made it possible to adopt a time-like iteration scheme- for solving the governing equations. The success of 
this approach prompted a series of attempts to develop preconditioning methods to be used in conjunction with 
compressible Navier-Stokes equations (Turkel [32], Van Leer et al. [33], Choi and Merkle [34], Venkataswaran 
and Merkle [35]). This method is especially valuable for solving aerodynamic flow fields involving multiple 
scales (i.e., flow fields involving regions of very low speed in generally high-speed flow problems). This 
formulation continues to be utilized in simulating multi-phase flow such as cavitation in liquid flow [36]. 

Numerical Schemes 

It is common to see that the artificial compressibility code resembles that of a compressible flow. For example, to 
alleviate the difficulty of solving an unfactored implicit scheme, indirect methods have been devised by many 
researchers. The alternating direction implicit (ADI) scheme by Beam and Warming [37] or Briley and McDonald 
[38] approximates the implicit operator in the unfactored scheme by a product of three one-dimensional operators. 
This scheme is unconditionally stable in two dimensions. In three dimensions, with the use of numerical 
dissipation, this scheme becomes conditionally stable. The ADI scheme introduces a factorization error that 
reduces the rate of convergence. This method introduced perhaps the first implicit scheme to CFD, thus viscous 
flow computation became more feasible. With the help of reduced cost by the diagonalization of the Jacobian 
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(Pulliam and Chaussee [39]), this scheme has been used in the development of many flow solvers. At NASA 
Ames, INS3D, the first incompressible Navier-Stokes solver in a generalized three-dimensional coordinate, 
(Kwak et al. [27]), was developed using this scheme. The original INS3D code was extensively used as a primary 
design tool for an upgraded design of the Space Shuttle main engine (SSME) hot gas manifold by Chang et al. 
[28]. This was probably the first major CFD contribution to rocket propulsion system design. 

In the artificial compressible formulation, since the governing equations are changed into the hyperbolic- 
parabolic type, some of the upwind differencing schemes developed for the compressible Euler and Navier- 
Stokes equations [40-43] can be used. These schemes were implemented into the flow solver successfully using 
the artificial compressibility method by Rogers et al. [3 1]. A variety of iterative schemes can be implemented in 
conjunction with the upwind schemes. For example, Rogers et al. [31] solved the resulting set of numerical 
equations using an unfactored line relaxation scheme similar to that employed by MacCormack [44]. 

2.4 Other Methods 

Although the evolution of incompressible flow solver development has been explained here using examples the 
authors of this paper are directly involved with, a large number of research projects have been ongoing in the 
incompressible CFD community in general. In particular, development along the finite element or unstructured 
grid approach has not been examined in this paper. Readers interested in this area are referred to a fairly recent 
book by Gresho and Sani [9], and a number of articles in the book by Hafez [10] (for example, see articles by 
Glowinski [45], Morgan [46], Tezduyar [47], as well as a many articles covering numerical methods and special 
topics). Several other methods related to incompressible flow computations can be found in the references cited. 

In general it is the authors’ opinion that there are no longer any major differences in the basic methodologies 
and tools used by different organizations. Differences in robustness, computational efficiencies, and accuracy 
usually come from how the solution procedures are set up in conjunction with geometry and grid handling, 
software engineering, and some subtle algorithmic aspects not well-suited to presentation in written papers. 

3. Computed Examples 

Solutions to some basic problems 

Validation computations of solution algorithms and flow solvers have been reported in the literature over the 
years; e.g., problems like square cavity, flow over a circular cylinder, duct flow of various shapes, backward 
facing steps, juncture flow, and tip and wake vortex flow. Even though the geometry of these flow problems are 
relatively simple, they offer rich sources for studying flow physics and opportunities to prove the validity of 
particular numerical methods and procedures. One simple example is shown in Figure 1, where the artificial 
compressibility and pressure projection methods are compared in order to evaluate the capability for predicting 
unsteady phenomena. Once the final solutions are obtained, the two approaches can produce equivalent results, 
however, the iterative processes for the artificial compressibility method and time advancing scheme for the 
pressure projection scheme have to be properly handled [48]. 
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Figure 1. Evolution of the stagnation point for an impulsively started flat plate: Comparison of artificial 
compressibility (denoted by INS 3D -UP), projection method (INS3D-FS), finite-element solution [49], and water tank 
experiment [50], 

Validation for real-world applications 

Validation of the solution procedure for real-world applications takes another' level of effort because the flow field 
involving complex geometry often requires high resolution locally to resolve the complicated flow, as well as 
requiring evaluation of physical model being implemented. Such a validation effort is of significant importance to 
problems like submarines with propulsor and for turbopump subsystems for liquid rocket engines. 

Figure 2 illustrates one such validation using a pump impeller. This validation is the first step toward further 
application of the simulation procedure to an entire turbopump system. Figure 2a shows an impeller geometry 
with surface grids. This particular configuration was designed to collect experimental data for developing an 
advanced impeller. Figure 2b shows some comparison of computed results with test data. Even with this 
simplified geometry, testing, which involves variations in flow conditions and shroud geometry, is not 
straightforward. Simulation involves geometric modeling and grid generation, as well as approximations related to 
numerical methods and physical modeling. As shown in the figures, the computed results are not on top of test 
data — -unlike the basic flow problems cited above. In engineering applications, especially for preliminary design 
applications, one needs to determine the acceptable level of accuracy. More importantly, the numerical simulation 
procedure must produce consistent results when flow conditions are varied as well as for variations in geometry 
such as different angle, length, and lean of impeller blades. 



Figure 2a. Geometry of a test impeller. 


Averaged Cm vs. Refative X: R/Rtip = 1.0275 



Figure 2b. Comparison of circumferentially averaged meridional velocity at the impeller exit. 

Applications to large engineering problems 
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Once validated, the computational procedure can then be extended to solve engineering problems. One of the 
most significant applications of incompressible flow simulation began in the early 1980s, when NASA carried 
out a series of upgrades to the Space Shuttle main engine (SSME) developed in the 1970s. One such effort was 
the powerhead redesign. Considered the backbone of the engine, the powerhead consists of the main injector 
assembly and pre-burners. Partially burned hot gas passes through Hot Gas Manifold (HGM) to the main 
injector assembly. From 1983-84, the powerhead redesign was undertaken, focusing on a two-duct HGM. 
NASA Ames and Rocketdyne collaborated in applying the INS3D code to this task. Since the Mach number of 
the partially burned gas was veiy low, the incompressible approximation was a valid assumption. The team of 
researchers successfully applied a CFD simulation procedure to this task for enhancing the performance of the 
SSME powerhead. The new design was developed based on a two-duct hot gas manifold. The two-duct design 
replaced the previous three-duct engine, smoothing the fuel flow, reducing pressure and turbulence, and 
lowering temperatures in the engine during operation. This two-duct design first flew on the shuttle in July 
1995. It significantly improved fluid flow in the system, thus reducing maintenance and enhancing overall 
engine performance. This pioneering work was probably the first major application of incompressible CFD to a 
rocket propulsion system (see [28] for more detail). 


Since this work, many incompressible flow software packages have been developed for analyzing complex 
engineering systems. One such example is the UNCLE code [36], which has been applied to various naval 
hydrodynamic applications. Many other equally significant software packages exist, as well. It is not the 
purpose of this paper to list the inventory 7 of such codes. 


A typical process of flow simulation, especially, for high-fidelity unsteady flow, requires large amounts of both 
computing Lime and huitiaii tune in problem setting and data processing. A substantial reduction in 
computational time for 3-D unsteady flow simulations is desired to reduce the design-cycle time of, for 
example, a pump system. Part of this speedup will be due to enhancements in computer hardware. The 
remaining portion of the speedup must be contributed by advances in grid-generation procedures, flow solution 
algorithms, and by efficient parallel implementations. In the following example, an overset strategy [51] was 
adopted to generate the grid for multiple components of a pump. To reduce human and computing time for grid 
generation, scripting [52] has been developed that substantially reduce the grid-generation portion of the total 
simulation time. Figure 3 illustrates this procedure. 


INLET GUIDE VANES (N repeated blades, no tip clearance) 



Figure 3. Illustration of automated grid-generation procedure with turbopump scripts for inlet guide vanes. 

Parallel computing strategies vary depending on computer architecture such as memory arrangement relative to 
processing units. Two approaches are illustrated here: the first approach is hybrid MPI/OpenMP and the second 
is Multi-Level Parallelism (MLP), developed at NASA Ames. The foimer is obtained by using message-passing 
interface (MPI) for inter-zone parallelism, and by using OpenMP directives for intra-zone parallelism. INS3D- 
MPI is based on the explicit message-passing interface across MPI groups and is designed for coarse-grained 
parallelism. 
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Figure 4. Parallel performance of INS3D-MLP on SGI Origin platform. 

The second approach to parallel computing is obtained by using NAS-MLP routines. This approach differs 
from the MPI/OpenMP approach in a fundamental way in that it does not use messaging at all. All data 
communication at the coarsest and finest level of parallelism is accomplished via direct memory referencing 
instructions. This approach also provides a simpler mechanism for converting legacy code than does MPI. For 
shared memory MLP, the coarsest level parallelism is supplied by spawning independent processes via the 
standard UNIX fork. The boundary data for the overset grid system is archived in the shared memory arena by 
each process. Other processes access the data from the arena as needed. Figure 4 shows the MLP scalability 
for impeller computations using 19.2 million grid points. 

Biomedical / Biofluid Applications 

Extension of incompressible flow methods to blood flow has been of interest to biomedical researchers for 
many years. For the cardiovascular system, the brain, and other parts of the human body, branching of blood 
vessels involves bifurcations, most of which are non-symmetric. Therefore, bifurcation has been a popular 
subject for blood flow simulations. Bifurcation problems offer the opportunity to address numerical issues 
involving grid generation, as well as to study basic fluid dynamic phenomena relevant to blood circulation 
simulation. When the size of a blood vessel becomes very small (as in capillaries), non-Newtonian 
characteristics become significant, thus requiring physical modeling. In addition, the vessel wall is elastic, 
requiring a structural model to account for geometric changes, depending on wall stresses. Early pioneering 
work on this subject was done by Womersley [53]. More sophisticated computations have since been performed 
for regular and non-regular bifurcating arteries including stenotic vessels (see Berger [54] for a review). 

The human circulatory system is like a huge tree with many branches of various sizes. Therefore, many 
computational studies have been performed using a truncated geometric model. One difficulty of simulating a 
truncated arterial system results from setting proper boundary conditions, especially at the downstream 
boundary. To account for the large arterial network, Quarteroni and his colleagues have developed a circuit 
analog, which has been applied to more inclusive circulatory systems modeling such as the cardiovascular 
system and the Circle of Willis in the brain (see [55] among many publications). An alternative way of 
imposing downstream boundary conditions has been developed by Tim David and his colleagues [56] in their 
computation of the Circle of Willis (CoW). 

One crucial step in biomedical flow simulation is to construct computational geometry starting from anatomical 
data. Kim et al. reconstructed a 3-D, anatomically realistic CoW geometry from human- specific magnetic 
resonance angiography (MRA) provided by researchers at Canterbury University in New Zealand (Tim David, 
Mechanical Engineering Dept; John Fink and Mike Hurrell, Christchurch Hospital). The image segmentation 
procedure is illustrated in Figure 5. The raw MR images were converted to the RGB graphic file format. After 
extracting the segments of interest by filtering voxels with intensities below a certain threshold, a segment- 
outlining algorithm was used to display the extracted objects on each sectional layer. 
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Figure 5. Image segmentation from a magnetic 
resonance image for a human-specific Circle of Willis. 
(MRA provided by John Fink and Mike Huirell, 
Christchurch Hospital, New Zealand.) 



Figure 6. Overlapping of a pulsatile blood flow 
simulation at peak systole on a magnetic 
resonance angiogram. 


With a non-Newtonian blood model, a model for deformable blood vessels, truncated downstream boundary 
conditions, and an auto-regulation model, Kim et aL [57] simulated unsteady blood circulation through this 
CoW under various gravitational conditions (a snapshot is shown in Figure 6, above). 


Modeling and simulation of the cardiovascular system has been the subject of many investigations because of 
its importance in biomedical research. Peter Hunter and his colleagues have been modeling cardiovascular 
systems, including multidisciplinary aspects, producing one of the most impressive results to date (e.g., see 
[58]). Earlier and perhaps the most elaborate physical model of the heart was pioneered by Peskin and his 
colleagues (see [59] for example) since 1970s. This model included blood, wall structures, and electrical field 
activating heart muscles. 


Another interesting application of blood flow simulation is related to artificial devices such as artificial hearts, 
ventricular assist devices (VADs). and heart valves. Because the demand for transplant organs far exceeds the 
number of donors, the need for artificial devices — to be used either as a temporary device or as a permanent 
replacement for a natural organ — becomes increasingly high. Accurate quantification of blood flow plays a 
crucial role in developing these devices. Thus, CFD simulation of blood flow in and around these artificial 
devices has become an indispensable part of the design. One such example is the recent CFD-aided design of 
the DeBakey VAD, which enabled human implantation by removing thrombus formation and lowering 
hemolysis to an acceptable level for human application. This effort set a new milestone for CFD application to 
the biomedical area [60]. 


All such blood flow computations may be regarded as a branch of incompressible CFD. The number of CFD 
applications to blood flow and biomedical problems are increasing rapidly, and the work cited here represents 
only small samples of the vast amount of ongoing and existing work. 

4. Challenges 

Incompressible CFD methodologies have been advanced along with computational technologies in general. 
Many fluid engineering problems can now be simulated, however, mostly at a single-component level. For 
example, it is possible to generate solutions to problems like a turbopump (inducer-impeller-diffuser), a naval 
vehicle at a steady motion including propulsor by a model, and a truncated model of the brain or heart. To 
realize the full benefits of CFD, more inclusive modeling will be required, such as, systems of pumps including 
multiples of pumps and feed lines, vehicle in maneuver with propulsor, a complete or more inclusive human 
circulatory system from the heart through the aorta and all the way to the brain and kidney capillaries. Attempts 
to solve these types of problems have been made with some qualitative successes. However, the predictive 
capability is still very limited and prediction with accurate physics is yet be accomplished. This will require 
inclusion of not only fluid dynamic modeling but modeling of other quantities like thermal loading, structural 
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properties (such as structural behavior of arterial walls), and cavitation physics. These computations will require 
not only large computing resources but large data storage and management technologies, as well. 

Flow solver codes and software tools have been developed to the point that many daily fluid engineering 
problems can now be computed routinely. Some of the physical models, such as those for turbulence and 
transition, however, have not advanced much since the ’70s or ’80s. Other models, like the cavitation model, 
have yet to be advanced to produce quantitative results for engineering. In aerospace design, the most 
productive aspect of CFD applications has been to predict relative change among design variations. To push the 
limit of operation and to try bold new ideas, more predictive capabilities will be needed for complicated flows 
involving transient phenomena, separation, tip vortex, and cavitation. For example, without accurate prediction 
capabilities for such quantities as cavitation and damaging frequencies, back flow, and rotational stall, CFD 
can’t be of much help in the development of an advance turbopump system. To make these advances, high- 
fidelity computation using high-end computing facilities are still be a “must” — despite the current euphoria 
about the PC clusters and grid computing. 

One final note is related to human resources. Even though CFD has advanced remarkably, many challenging 
cases require CFD experts. Computer science can automate a good portion of the CFD simulation processes, 
thus saving much human time required to obtain solutions. However, blind application of tools without 
understanding capabilities and limitations of the methods involved could lead to catastrophic engineering 
results. As in many other engineering and science disciplines, CFD researchers and practitioners need to 
understand the physics and engineering systems being simulated. Future experts will need to be cultivated who 
are willing to think through the flow physics in addition to the software engineering. 
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